Motion-induced synchronization in metapopulations of mobile agents 
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We study the influence of motion on the emergence of synchronization in a metapopulation of 
random walkers moving on a heterogeneous network and subject to Kuramoto interactions at the 
network nodes. We discover a novel mechanism of transition to macroscopic dynamical order induced 
by the walkers' motion. Furthermore, we observe two different microscopic paths to synchronization: 
depending on the rules of the motion, either low-degree nodes or the hubs drive the whole system 
towards synchronization. We provide analytical arguments to understand these results. 
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The spontaneous emergence of synchronization in sys- 
tems of coupled dynamical units [l|, Q underlies the de- 
velopment of coordinated tasks as diverse as metabolic 
cycles in eukaryote cells, cognitive processes in the hu- 
man brain and opinion formation in social systems 0-Q ■ 
In the last decade complex networks theory has revealed 
that the topology of the interactions in a complex sys- 
tem has important effects on its collective behavior [3, 111]. 
As a consequence, many recent studies have considered 
dynamical systems coupled through non-trivial topolo- 
gies @, uncovering the impact of the structure of the 
network on the existence l(j- 14 1 and stability [H - 18 1 of 
synchronized states. Quite frequently, the interactions 
among the units of a complex system keep changing over 
time. Their evolution can either be driven by the syn- 
chronization process itself, as in models of coevolving net- 
wor ks [mil], or be determined by the fact that each 
unit moves at random over a continuous and homoge- 
neous space and interacts only with other units within a 
given distance jioT - Eril ]. In many cases, the motion of the 
agents takes place on discrete and heterogeneous media, 
that can be represented as complex networks. Typical 
examples include users browsing the World Wide Web, 
airplane passengers traveling throughout a country, or 
people playing online social games |27H29|. In such sys- 
tems, both the rule of motion adopted by the agents, and 
the heterogeneity of the environment, have an impact on 
the emergence and stability of collective behaviors. For 
this reason, metapopulation modeling has been success- 
fully employed to explore the combined effect of mobility 
and non-trivial interaction patterns in different contexts, 
including the study of epidemic spreading and chemical 
reactions [13, [H, Wt&fy . 



In this Letter we propose a metapopulation model to 
study the emergence of synchronization in populations 
of individuals moving over discrete heterogeneous envi- 
ronments, and interacting through nonlinear dynamical 



equations. We assume that each agent is characterized 
by an internal state (or opinion), and that it moves over 
the environment trying to synchronize its opinion with 
that of the other individuals. Thus, the evolution of the 
system is driven by the interplay of two concurrent pro- 
cesses: on one hand, the interaction of neighboring agents 
drives their internal state towards local consensus; on the 
other hand, agents' motion dynamically changes the pat- 
tern of interaction and allows each agent to be exposed to 
different opinions. We discover a novel mechanism of syn- 
chronization that we name motion-induced synchroniza- 
tion, since the transition from disorder to macroscopic 
order is controlled by the value of the parameter tuning 
the motion of agents. Furthermore, we show that there 
are two different microscopic mechanisms driving the sys- 
tem towards synchronization, according to whether the 
walkers prefer to visit or to avoid high-degree nodes. 

Our metapopulation model consists of two layers. At 
the bottom layer we have a set of W of mobile agents 
(walkers). Each agent i (a = 1,2,... W) is a dynamical 
system whose internal state at time t is described by a 
phase variable, 9i(t) € [0,27r)%, and changes over time 
as a result of the interactions with other agents. The top 
layer consists of a complex network with N nodes and 
E edges, which represents the environment into which 
agents interact (nodes) and move (edges). The network 
is described by an adjacency matrix A, whose entry ojj 
is equal to 1 if nodes / and J are connected by an edge, 
and otherwise (here and in the following we indicate 
nodes of the graph in uppercase letters, and walkers in 
lowercase). At any given time, each agent is located in 
one of the nodes of the network. The agent interacts for 
a fixed time interval with other agents at the same node, 
trying to synchronize its phase with the others'. Then, it 
moves to one of the neighboring nodes, chosen according 
to a one-parameter motion rule. More precisely, assume 
that at time t we have i £ I, i.e. agent i is at node I. The 
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evolution of the phase ftj (t) of agent i is ruled by an all- 
to-all Kuramoto-likc interaction with the other walkers 
being on the same node I at time t 0, HH [lf| : 



i(t) 



UJi 



i(t)),VieI (1) 



where lji is the internal frequency of agent i and A is 
a control parameter accounting for the strength of the 
interaction among walkers. Notice that, when the phases 
of the agents evolve according to Eq. ([T]), the system is 
not driven by a single global mean-field (as in the classical 
all-to-all Kuramoto model). Instead, each oscillator i 
interacts with the local mean-field phase due to all the 
oscillators being at the same node as i. At regular time 
intervals of length A all the agents perform one step of a 
degree-biased random walk on the network. Namely, we 
assume that a walker at node / moves to a neighboring 
node J with a probabili ty p roportional to the degree kj 
of the destination node [23 l2Sl : 
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Here a is a tunable parameter which biases agents' mo- 
tion either towards low-degree nodes (a < 0) or towards 
hubs (a > 0). For a = 0, we recover the standard (un- 
biased) random walk. In summary, the metapopulation 
model has three control parameters: A regulating the in- 
teraction strength among walkers, a tuning the rule of 
their motion, and A fixing the ratio between the time 
scales of interaction and motion. The degree of synchro- 
nization of the whole metapopulation at time t is mea- 
sured by the global order parameter: 
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where r ~ if the phases of the agents are completely 
incoherent, while r = 1 when the system is fully synchro- 
nized. In order to quantify the degree of synchronization 
of a single node / we introduce the local order parameter: 
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where wi(t) is the number of agents at node I at time 
t. When the phases of the walkers at node I are fully 
synchronized, the local order parameter of the node is 
equal to 1, while in the case of complete local disorder 
we have ri(t) = 0. We can also quantify the average local 
synchronization of the network r\ oc (t) as the average of 
ri(t) over all nodes, i.e.: 
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FIG. 1. (color online) Phase diagram r(a, A) of the metapop- 
ulation model on uncorrelated SF networks with 7 = 3 (panel 
A) and on the US air-transportation network (panel B) . Both 
networks have N = 500 nodes. The dashed curve in the first 
panel is a lower-bound analytical prediction for the onset of 
synchronization in uncorrelated graphs obtained from Eq. §9§ 
for a > — 1, and from Eq. (|10[) for a < — 1. Panel C shows 
r(a) for A = 0.08, corresponding to the horizontal line in 
panel A. 



We notice that, having rj(t) ~ 1 V/, or equivalently, 
ri oc (t) ~ 1 is a necessary but not sufficient condition to 
attain global synchronization. In fact, in the limiting 
case in which there is no motion (A — > oo) and A is large 
enough, it is possible to have r/(i) ~ 1 Vi and, at the 
same time, r(t) ~ 0. See Appendix. 

We have simulated the metapopulation model on vari- 
ous synthetic networks and, as an example of a real com- 
plex network, on the US air-transportation system, which 
includes the flight connections, between the N = 500 
largest airports in the US. Thanks to its intrinsic nature 
as a backbone for human transportation, the US airports 
network has already been used to investigate reaction- 
diffusion dynamics in metapopulation models (2~7| . This 
network has a long-tailed degree distribution, exhibits 
degree-degree correlations and is relevant for the present 
study because opinion formation in real social systems 
is often mediated by information, communication and 
transportation networks having similar structural prop- 



erties. We have generated, for comparison, uncorrelated 
scale-free (SF) networks with N = 500 nodes and a 
power-law degree distribution P(k) ~ fc -7 with a tun- 
able value of the exponent 7 |36[. For the numerical 
simulations of the model, the initial phases of the os- 
cillators have been sampled uniformly in [0,2-7r), and 
their internal frequencies from a uniform distribution 
g(ui) = 1/2 V uj G [—1, 1]. We started from a stationary 
distribution of W = 5000 walkers over the networks [37( , 
and we integrated Eq. (p} for all the agents for a time 
Iq = mA, where m = 10 4 is the number of random 
walk steps performed. After this transient, we estimated 
global and local synchronization parameters (r, ri oc , and 
rj for I = 1 , . . . , N) by respectively averaging the values 
obtained from Eq. ([3]), (|4|) and (O over a time window 
of length T = 2mA. In Fig. [JJwe show the global order 
parameter r as a function of the coupling strength A, and 
of the walker bias a. The two phase diagrams have been 
obtained setting A = 0.05. Qualitatively similar results 
have been obtained for different values of A (see Ap- 
pendix). The SF network reported in Fig. [TJA has 7 = 3, 
although similar results have been obtained for other val- 
ues of the exponent in the range 2 < 7 < 3. As expected, 
by increasing A at a fixed value of a, i.e. keeping fixed 
the rules of motion, we observe a phase transition from 
the incoherent phase (r ~ 0, dark regions of the dia- 
grams) to a synchronized state (r^0, bright regions of 
the diagrams). However, the precise value for the onset 
of synchronization, namely the critical value A c for which 
the incoherent state becomes unstable, strongly depends 
on the motion bias a. In particular, we find that X c (a) 
is first increasing as function of a, and then decreasing. 
The function A c (a) reaches its maximum A" lax at a par- 
ticular value of a, namely at a* ~ — 1 for the SF net- 
work (we have checked that this is true for any value of 
7 in the interval (2,3]), and at a* ~ —0.5 for the air- 
transportation network. As a consequence of the shape 
of A c (a), for a wide range of values of the strength A such 
that A < A™ 11 *, we observe a novel mechanism of motion- 
induced synchronization. This means that we can fix the 
value of the interaction strength A, and we can control 
whether the system is in the incoherent phase (r ~ 0, 
dark regions) or in the synchronized state (r 7^ 0, bright 
regions) solely by changing the rule of motion. Let us 
focus on the case of the SF network shown in Fig. [TJA 
and suppose to keep A fixed at 0.08 (see the horizon- 
tal line in the figure). As shown in Fig. [T]C, we can 
start with a value of a within the incoherent phase, for 
instance a = a* = — 1, and consider the behavior of 
the system as we decrease the value of a. When a gets 
smaller than a particular critical value a cl (A), in this 
case a ci (0.08) ~ —2.0, wc observe a transition from the 
incoherent to the synchronized phase. Conversely, we 
can start at a = — 1 and get a synchronized state by in- 
creasing the motion bias parameter to values larger than 
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FIG. 2. (color online) Panel A and B: the average local or- 
der parameter rk of nodes of degree k as a function of k, for 
various values of A, and for two fixed values of the bias, re- 
spectively q = —2.0 (panel A) and a — —0.25 (panel B), 
corresponding to the two vertical lines in Fig. [T]A. Panel C 
and D: rk for A = 0.08 and various values of a (respectively, 
a < — 1 in panel C and a > — 1 in panel D) corresponding to 
the horizontal line in Fig. [T]A. 



ling the agents' motion can effectively produce dramatic 
changes in the macroscopic synchronization state of the 
system. 

We have found that the bias in the motion affects the 
onset of synchronization also at the microscopic level. To 
illustrate this result we look at the microscopic paths to 
synchronization [l2j as we increase A, by following the 
two vertical lines shown in Fig. [TJA. In particular, we 
have computed the value of the local order parameter 
for each of the nodes of the graph, as in Eq. ([4]), and 
we have grouped nodes by degree classes. We computed 
the average value of the local synchronization of nodes of 
degree k as: 



1 N 
— ^V /( 5(fc/,fc) 
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-0.2. A fine tuning of the rules control- 



where Nk = NP(fc) is the number of nodes of degree fc. 
We report in Fig. ®A and [5]B the quantity divided 
by the value r^(A ~ 0) obtained when A is close to 
as a function of k, and for different values of A. Panel 
A corresponds to a = —2.0 and panel B to a = —0.25. 
When a = —2.0 the nodes having small degree are the 
first ones to attain local synchronization as soon as A 
crosses the critical value A c (— 2.0) ~ 0.08; conversely, for 
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a = —0.25, the hubs are the nodes which synchronize 
first when A > A c (— 0.25) ~ 0.07. We thus observe two 
microscopic paths to synchronization: cither driven by 
low-dcgrcc nodes (a < a*), or by the hubs (a > a*). 
The two different synchronization mechanisms are also 
evident by following the horizontal line in Fig. [T]A, i.e. 
by plotting for a fixed value of A and different values 
of a, as shown in Fig. [2lC andElD. 

The effects of motion on synchronization can be ex- 
plained by analytical arguments in the case of networks 
without degree-degree correlations. In particular we de- 
rive, as follows, a lower-bound estimate for the critical 
strength A c as a function of a. The average number 
wj of biased random walkers at a node / of an undi- 
rected connected graph without degree-degree correla- 
tions reads (37. 38|: 
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where cj = X)j=i a u^j- F° r a given a, the value wj de- 
pends only on the connectivity kj of the node, so that all 
the nodes with the same degree will have the same aver- 
age number of walkers. Thus, in the following we indicate 
as Wk the number of agents on a node with degree k. We 
consider now the two limiting cases A — > and A — > oo. 
When A — > (fast-switching approximation) the agents 
on a node interact for an infinitesimal time interval be- 
fore moving to another node. In this limit we have a well- 
mixed population of oscillators that can be approximated 
as a single all-to-all Kuramoto model of W elements. 
Thus, the critical value of the coupling A, in this case 
reads A c = 2/[W ■ ir ■ g(0)] [I El], and does not depend 
on a. When A — >• oo (slow-switching approximation), i.e. 
when the walk is much slower than the Kuramoto dynam- 
ics, each node of the network is an all-to-all Kuramoto 
system independent from the others. For a fixed value of 
A, in some nodes the oscillators will reach local synchro- 
nization before eventually walking away, while in some 
other nodes they will not. In fact, the critical coupling 
strength of a node / is that of a set of wi all-to-all coupled 
Kuramoto oscillators: A c (7) = 2/[wiir ■ g(0)]. Hence, the 
critical coupling strength for the local synchronization of 
the walkers at a node of degree k reads: 



Ac(fc) 
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where we have made use of Eq. ([7]). Therefore, in the 
slow-switching approximation, at fixed values of W/N, 
a and A, only agents at nodes of degree k such that 
A c (fc) < A will attain local synchronization. Conse- 
quently, a necessary but not sufficient condition to have 
global synchronization is that there is at least one node 
J in the graph for which X c (kj) < A. Eq. ([8]) sheds 
light on the two different microscopic paths to synchro- 
nization observed in Fig. [21 Let us indicate as fc m i n and 



fcmax respectively the minimum and the maximum de- 
gree in the network. Consider two values of a, one larger 
and one smaller than a* = —1, for instance the two 
values a = —2 and a = —0.25 corresponding to the 
two vertical lines in Fig. [T]A. If we start increasing A 
from A = 0, the slow-switching approximation predicts 
no local synchronization until A becomes larger than the 
smallest value of A c (fc), corresponding to k = fc m in if 
a < —1, or to k = fe max if a > — 1. At this point, if 
a < — 1 (a > —1) the walkers at nodes with the small- 
est (largest) degree attain local synchronization. If we 
keep increasing A, local synchronization is progressively 
reached also at nodes with larger (resp. smaller) degrees 
when a < — 1 (resp. a < —1). We can therefore de- 
rive a lower-bound A c (a) for the curve A c (a) delimiting 
the synchronization region in Fig. [T]A, by considering 
the smallest value of A at which at least one class of 
nodes attains local synchronization. In particular, for a 
SF network with P(k) ~ fc~ 7 , as the ones used in our 
simulations, we get: 
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when a > — 1, and: 
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for a < —1. We notice that if the graph is uncorrelated, a 
motion rule with a = — 1, leads to a uniform distribution 
of the walkers over the nodes, since wi = W/JWZ in 
Eq. ([7]), and all the nodes attain local synchronization 
altogether at A = 4N/ (Wir), as can be seen from Eq. ([8]). 
This corresponds to the largest possible value A™ ax of 
the critical interaction strength. In Fig. [T]A wc report 
as dashed line the curve A c (a) obtained for the same 
values of W/N, fc m i n , W and 7 used in the numerical 
simulations. Although the slow-switching approximation 
provides only a lower-bound for the critical interaction 
strength, it works quite well for networks without degree- 
degree correlations, and it also predicts quite accurately 
the position of the cusp at a = a* = — 1 for any value of 
7 in (2,3]. The good agreement between the prediction 
in the slow-switching approximation and the numerical 
simulations indicates that a value A = 0.05 corresponds 
indeed to an intermediate regime of the system. In this 
regime the motion is fast enough to allow a certain level 
of mixing and the attainment of global synchronization 
and, at the same time, it is slow enough to avoid full 
mixing, for which A c (a) = A c Va (see Appendix). 
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Appendix 

In our model, the parameter A controls how often the 
agents perform a step of random walk. By tuning A one 
can change the ratio between the time-scale of the in- 
teraction and that of the walk. Here we study how the 
value of A affects the synchronization of the system. We 
first consider the model in the limiting case A — > oo, 
in which the agents are not allowed to move. For each 
value of a, we distributed a population of W = 5000 
walkers across the nodes of a random scale-free network 
with N = 500 nodes and P(k) ~ fc -3 , according to the 
stationary distribution of Eq. (7). Since there is no mo- 
tion, each oscillator will remain at the initial node and 
will interact with the same set of oscillators for all the 
duration of the simulation. In this limit, the metapop- 
ulation model is equivalent to a set of ./V independent 
all-to-all Kuramoto systems, with the system at a node 
of degree k having a critical interaction strength given 
in Eq. (8). In Fig. IS- II we report the local and global 
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FIG. S-2. Cross-section plots of the two phase diagrams in 
Fig. IS- 1 1 We report r\ oc and r as a function of the coupling 
strength A, for three different values of a. File absence of 
motion hinders global synchronization (r ~ 0, orange dashed 
lines), even if, for an appropriately large value of A, all the 
nodes of the network can reach complete local synchronization 
(floe = 1, solid lines). 
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FIG. S-l. Phase diagram reporting the local and global or- 
der parameters, ri oc (upper panel) and r (lower panel), as 
a function of a and A for the metapopulation model in the 
limit A — > oo (W = 5000 agents on a scale-free network with 
TV = 500 nodes). When there is no motion, the system is 
globally incoherent, even if the local synchronization at the 
nodes can be enhanced at will by increasing the value of the 
interaction strength A. 



order parameters, n oc and r, as a function of a and A. 
We observe that, for each value of a, there exists a criti- 
cal value of A such that at least one node of the network 
can attain local synchronization, and by increasing A we 
can reach high values of r\ oc (panel A). Conversely, the 
global order parameter r always remains close to zero, 
and no global synchronized state is found for any value 
of a and A (panel B). Notice that the phase diagram of 
Fig. IS-1I B looks quite different from the one reported in 
Fig.lA of the Letter, which corresponds to a simulation 
with A = 0.05 on the same network. This indicates that, 
in the absence of motion, the system will remain incoher- 
ent at a global scale, even if synchronization can emerge 
at the level of network nodes. The behavior of the model 
for A — > oo is better illustrated by the cross-section plots 
shown in Fig. IS-21 Here, we report the values of local 
and global order parameters, ri oc and r, for three dif- 
ferent choices of the bias, namely for a = —1.5 (black), 
a = —1.0 (red) and a = —0.5 (blue). Notice that if A 
is large enough all the nodes will eventually attain local 
synchronization (rioc — l)j but for any combination of a 
and A, the system remains globally incoherent (r ~ 0). 
Notice that, as expected, all the nodes achieve full local 
synchronization altogether when a = —1.0, i.e. when the 
system is initialized with an equal number of walkers at 
each node. 

We now consider the metapopulation model for finite 
values of A. In Fig. lS-3l we report the value of the global 
order parameter rasa function of A, for A = 0.1 and two 
values of the motion bias, namely a = —2.0 (red open cir- 
cles) and a = 0.5 (blue filled circles). When the value of 
A is large, i.e. when the motion is rare with respect to 
agents interaction, the global order parameter decreases 



dramatically and approaches the behaviour of the lim- 
ing case A — > oo. This means that even if the value of 
A is large enough to guarantee that all the agents on a 
node will attain full synchronizaton between two subse- 
quent steps of the random walk, the poor mixing due 
to rare motion prevents the emergence of global order. 
Conversely, if A is small then the interaction interval at 
each node is not large enough to allow local synchroniza- 
tion; nevertheless, the presence of fast motion enhances 
mixing and promotes the convergence of each oscillator 
towards a global synchronized state. 




FIG. S-3. Global order parameter r as a function of A for 
A = 0.1 and two values of the motion bias, respectively 
a — —2.0 (red open circles) and a — 0.5 (blue filled cir- 
cles). Global synchronization is attained when agents move 
frequently (A — > 0), while the system remains in an incoher- 
ent state if the motion is too slow (A — > oo). 



